2023 Dev/GlowwormScaleRanking.R

In Hannahglover/Glowworm: Mapping genes to single cell RNAseq

GlowwormCompRanking = function(Input, ShowTop = "All", Colors = "Default", RotateXText = 0, LabelSize = 11, Legend=T, Frequency = T, Flip = F){
  
  if('Default' %in% Colors & "Dataset" %in% colnames(Input@MetaData) ==T & length(unique(Input@MetaData[["Dataset"]])) > 1){Cols = hue_pal()(length(unique(Input@MetaData[["Dataset"]])))
  }else if('Default' %in% Colors & "Dataset" %in% colnames(Input@MetaData) == T & length(unique(Input@MetaData[["Dataset"]])) == 1 | 'Default' %in% Colors & "Dataset" %in% colnames(Input@MetaData) == F){Cols = "#808080"
  }else if(class(Colors) == "character" & length(Colors) < length(unique(Input@MetaData[["Dataset"]]))){stop("Need to provide the same number of colors as datasets , or use 'Default'")
  }else if(class(Colors) == "character" & length(Colors) >= length(unique(Input@MetaData[["Dataset"]]))){
    Cols = Colors
  }
  if(RotateXText == 90){
    SetVjust = 0.5
    SetHjust = 1
  }else if(RotateXText == 45){
    SetVjust = 1
    SetHjust = 1
  }else{
    SetVjust = 0
    SetHjust =  0.5
  }
  
  
  #CalculateDeltas = as.data.frame(rowMeans(Input@GlowwormScaleOutput))
  #colnames(CalculateDeltas) = "TargetGenes"
  #CalculateDeltas$Secondary = rowMeans(Input@GlowwormScaleOutput_Secondary)
  #CalculateDeltas$Delta = CalculateDeltas$TargetGenes-CalculateDeltas$Secondary
  #CalculateDeltas$Populations =  row.names(CalculateDeltas)
  #CalculateDeltas$Populations =  gsub("\\^", " - ", CalculateDeltas$Populations)  
  
  #Input@GlowwormScaleAllCells$Combined = NULL
  #Input@GlowwormScaleAllCells_Secondary$Combined = NULL
  #  AllData = Input@GlowwormScaleAllCells
  #  AllData$Combined = NULL
  
  #  CalculateDeltas = as.data.frame(rowMeans(AllData))
  #  colnames(CalculateDeltas) = "TargetGenes"
  #  CalculateDeltas$Secondary = rowMeans(Input@GlowwormScaleAllCells_Secondary)
  #  CalculateDeltas$Delta = CalculateDeltas$TargetGenes-CalculateDeltas$Secondary
  #  CalculateDeltas2 = merge(CalculateDeltas, Input@MetaData, by = 0)
  #  CalculateDeltas2$Combined =  paste(CalculateDeltas2[["Subclass"]], CalculateDeltas2[["Class"]], sep="^")
  #  CalculateDeltas2$nFeature_RNA = NULL; CalculateDeltas2$nCount_RNA = NULL; CalculateDeltas2[["Class"]] = NULL; #CalculateDeltas2[["Subclass"]] = NULL
  
  CompileTRes = as.data.frame(matrix(ncol = 3, nrow=0))
  colnames(CompileTRes) = c("Combined", "P.val", "meanDiff")
  i = 0
  AllPops = round(length(unique(Input@GlowwormScaleAllCells$Combined))/40)
  Progress =seq(1,length(unique(Input@GlowwormScaleAllCells$Combined)), AllPops)
  for(Comb in unique(Input@GlowwormScaleAllCells$Combined)){
    GetTargets = subset(Input@GlowwormScaleAllCells, Input@GlowwormScaleAllCells$Combined == Comb)  
    GetTargets2 = GetTargets %>% group_by(Barcs) %>% dplyr::summarise("AvgbyCell" = mean(Data))
    GetSecondary = subset(Input@GlowwormScaleAllCells_Secondary, Input@GlowwormScaleAllCells_Secondary$Combined == Comb)  
    GetSecondary2 = GetSecondary %>% group_by(Barcs) %>% dplyr::summarise("AvgbyCell" = mean(Data))
    Tres = t.test(GetTargets2$AvgbyCell, GetSecondary2$AvgbyCell, paired = TRUE, alternative = "two.sided")
    TresCompile = as.data.frame(t(c(Comb, as.numeric(Tres[[3]]), as.numeric(Tres[[5]]))))
    colnames(TresCompile) = c("Combined", "P.val", "meanDiff")
    CompileTRes = rbind(CompileTRes, TresCompile)
    i = i+1
    if(i %in% Progress){cat("=")}
  }
  CompileTRes$P.val = as.numeric(CompileTRes$P.val)  
  CompileTRes$meanDiff = as.numeric(CompileTRes$meanDiff)  
  
  # GetTargets3 = Input@GlowwormScaleAllCells %>% group_by(Combined) %>% dplyr::summarise("AvgbySubclass" = mean(Data))
  
  #if("Dataset" %in% colnames(Input@MetaData)){
  #}
  
  GeneData_t_Binary = as.data.frame(ifelse(Input@GlowwormScaleOutput > 0, 1, 0))
  GeneData_t_Binary[is.na(GeneData_t_Binary)] = 0
  SumBinary = as.data.frame(rowSums(GeneData_t_Binary))
  colnames(SumBinary) = "Frequency"
  
  SumBinary$GenesExpressed=apply(GeneData_t_Binary,1,function(x) paste(names(GeneData_t_Binary)[which(x==1)], collapse="/"))
  Output = list()
  for(y in c("Expression", "Add", "Product", "DivScaled", "Div")){
    CellNorm = y
    TopDeltas = merge(CompileTRes, SumBinary,by.x = "Combined", by.y = 0)
    
    if("Dataset" %in% colnames(Input@MetaData) == T){
      MetaData = Input@MetaData  
      MetaData$Combined = paste(MetaData$Subclass, MetaData$Class, sep="^")
      MetaData.red = MetaData %>% group_by(Combined, Dataset) %>% dplyr::summarise(n = n(), .groups = 'drop') 
      TopDeltas = merge(TopDeltas, MetaData.red, by = "Combined")  
    }
    TopDeltas$Significant = ifelse(TopDeltas$P.val < 0.05 & TopDeltas$meanDiff > 0.2, "Enriched", "")
    TopDeltas$meanDiff_Scaled = (TopDeltas$meanDiff- min(TopDeltas$meanDiff))/(max(TopDeltas$meanDiff) - min(TopDeltas$meanDiff))
    TopDeltas$Freq_Scaled = (TopDeltas$Frequency- min(TopDeltas$Frequency))/(max(TopDeltas$Frequency) - min(TopDeltas$Frequency))
    if(CellNorm == "Expression"){
      TopDeltas = data.frame(TopDeltas[order(-TopDeltas$meanDiff),])
      TopDeltas$RankScore = TopDeltas$meanDiff
    }else if(CellNorm == "Add"){
      TopDeltas$RankScore =   TopDeltas$Freq_Scaled + TopDeltas$meanDiff_Scaled
      TopDeltas = data.frame(TopDeltas[order(-TopDeltas$RankScore),])
    }else if(CellNorm == "Product"){
      TopDeltas$RankScore =   TopDeltas$Freq_Scaled * TopDeltas$meanDiff_Scaled
      TopDeltas = data.frame(TopDeltas[order(-TopDeltas$RankScore),])   
    }else if(CellNorm == "DivScaled"){
      TopDeltas$RankScore =   TopDeltas$Freq_Scaled/TopDeltas$meanDiff_Scaled
      TopDeltas = data.frame(TopDeltas[order(-TopDeltas$RankScore),]) 
    }else if(CellNorm == "Div"){
      TopDeltas$RankScore = TopDeltas$Frequency/TopDeltas$meanDiff
      TopDeltas = data.frame(TopDeltas[order(-TopDeltas$RankScore),])    
    }
    
    TopDeltas$Freq_Scaled = NULL; TopDeltas$meanDiff_Scaled = NULL
    
    TopDeltas2 = subset(TopDeltas, TopDeltas$Significant == "Enriched")
    
    if(ShowTop == "All"){
      TopDeltas2 = TopDeltas
      TopDeltas2$Significant = factor(TopDeltas2$Significant, levels = c("Enriched", ""))
      TopDeltas2 = data.frame(TopDeltas2[order(TopDeltas2$Significant, -TopDeltas2$RankScore),])    
    }else if(class(ShowTop) == "numeric"){
      TopDeltas2 = head(TopDeltas2, ShowTop)
    }else if(class(ShowTop) == "character"){
      MakeNum = as.numeric(ShowTop)
      if(class(MakeNum)== "numeric"){
        TopDeltas2 = head(TopDeltas2, MakeNum)
      }}else if(ShowTop == "Significant"){
        TopGenes2 = dim(TopDeltas2)[1]
      }else{stop("ShowTop must be 'All', 'Significant' or the number of populations you wish to plot ")}
    
    
    if(dim(TopDeltas2)[1] == 0){
      cat("\nNo populations were statistically enriched for the gene list\n")
    }else if(dim(TopDeltas2)[1] == 1){
      cat("1 population was statistically enriched for the gene list\n")
    }else if(dim(TopDeltas2)[1] > 1){  
      cat(paste("\n", dim(TopDeltas2)[1], " populations were statistically enriched for the gene list\n", sep=""))
    }
    
    if(Flip == T){
      TopDeltas2$Combined= factor(TopDeltas2$Combined, levels = rev(unique(TopDeltas2$Combined)))
    }else if(Flip == F){
      TopDeltas2$Combined= factor(TopDeltas2$Combined, levels = unique(TopDeltas2$Combined))
    }
    
    if(length(Cols) == 1){ 
      if(Frequency == T){
        FreqPlot = ggplot(TopDeltas2, aes(x=Combined, y=RankScore, size = Frequency))
      }else if(Frequency == F){
        FreqPlot = ggplot(TopDeltas2, aes(x=Combined, y=RankScore)) #+ theme_classic()  + geom_point( alpha = 0.7, position=position_jitter(w = 0.05, h = 0.02)) + theme( axis.text.x = element_text(angle = 90, hjust = 1, vjust=0.5, size=LabelSize), axis.text.y = element_text(size=LabelSize), legend.text = element_text(size=LabelSize), axis.title.y = element_text(size=LabelSize),  legend.title = element_blank()) +xlab("")+ylab("Cumulative Normalized Counts") +scale_size(limits = c(0, max(TopDeltas2$Frequency)))  + scale_color_manual(values = Cols)
      }
      
    }else if(length(Cols) > 1){ 
      if(Frequency == T){
        FreqPlot = ggplot(TopDeltas2, aes(x=Combined, y=RankScore, size = Frequency, color=Dataset)) 
      }else if(Frequency == F){
        FreqPlot = ggplot(TopDeltas2, aes(x=Combined, y=RankScore, color=Dataset)) 
      }
    }
    FreqPlot = FreqPlot + theme_classic()  + geom_point( alpha = 0.7, position=position_jitter(w = 0.05, h = 0.02)) + theme( axis.text.x = element_text(angle = 90, hjust = 1, vjust=0.5, size=LabelSize), axis.text.y = element_text(size=LabelSize), legend.text = element_text(size=LabelSize), axis.title.y = element_text(size=LabelSize),  legend.title = element_blank()) +xlab("")+ylab("Cumulative Normalized Counts") +scale_size(limits = c(0, max(TopDeltas2$Frequency)))  + scale_color_manual(values = Cols)
    
    if(Legend == F){
      FreqPlot = FreqPlot + NoLegend()    
    }
    if(Flip == T){
      FreqPlot = FreqPlot + coord_flip()
    }    
    
    Output[[y]] = FreqPlot
  }  
  #setOldClass(c("gg", "ggplot"))
  #setClass("GlowwormObj", representation(Plot = "ggplot", Data = "data.frame"))
  #Output = new("GlowwormObj", Plot = FreqPlot, Data = TopDeltas)  
  
  return(Output)
  
}